Tarea 1: Foundations

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Elaboración de Código

Objetivo

a la primera tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en análisis exploratorio de datos y conceptos introductorios de probabilidades. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de ejercicios prácticos con el fín de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Videos de las clases:

Elaboración de código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Pregunta 1: Visualización de Datos

Para esta pregunta usted deberá trabajar en base al conjunto de datos hearth_database.csv, el cual esta compuesto por las siguientes variables:

En base al dataset propuesto realice un análisis exploratorio de los datos (EDA). Para su análisis enfoquen el desarrollo en las siguientes tareas:

Respuesta

data <- read.csv("hearth_database.csv")

#summary(data

for (i in 7:14) {
  cat(colnames(data)[i],"\n")
  cat("Promedio: ", mean(data[,i]), "\n")
  cat("Mediana: ",median(data[,i]), "\n" )
  cat("Cuantiles: ", quantile(data[,i]), "\n")
  cat("Máximo: ", max(data[,i]), "\n")
  cat("\n")
}
cor(data[,7:14])
boxplot(x=data[,7:14],main="Boxplot data")
data <- read.csv("hearth_database.csv")
data$target <- as.numeric(data$target == "YES")

hist(data[,7][data$target == 1],nclass=10, main = "slope")
hist(data[,7][data$target == 0],nclass=10, main = "slope")

hist(data[,8][data$target == 1],nclass=10, main = "ca")
hist(data[,8][data$target == 0],nclass=10, main = "ca")

hist(data[,9][data$target == 1],nclass=10, main = "thal")
hist(data[,9][data$target == 0],nclass=10, main = "thal")

hist(data[,10][data$target == 1],nclass=10, main = "age")
hist(data[,10][data$target == 0],nclass=10, main = "age")

hist(data[,11][data$target == 1],nclass=10, main = "trestbps")
hist(data[,11][data$target == 0],nclass=10, main = "trestbps")

hist(data[,12][data$target == 1],nclass=10, main = "chol")
hist(data[,12][data$target == 0],nclass=10, main = "chol")

hist(data[,13][data$target == 1],nclass=10, main = "thalach")
hist(data[,13][data$target == 0],nclass=10, main = "thalach")

hist(data[,14][data$target == 1],nclass=10, main = "oldpeak")
hist(data[,14][data$target == 0],nclass=10, main = "oldpeak")

Pregunta 2: Teorema Central del Limite

Pruebe el teorema central del limite aplicando un muestreo de la media en las distribuciones Gamma, Normal y una a su elección. Grafique los resultados obtenidos y señale aproximadamente el numero de muestreos necesarios para obtener el resultado esperado, pruebe esto con las siguientes cantidades de muestreo \(\{10,100,1000,5000\}\). ¿El efecto ocurre con el mismo número de muestreo para todas las distribuciones?.

Respuesta

# Definición de variables o estructuras necesarias para el muestreo.
n <- 1000
tamaño <- 10
# Normal
sigma <- sqrt(5)
mu <-5
# Gamma
alfa <- 2
beta <- 1.5
# Poisson
lambda <- 5
#se inicializan vectores donde seran guardados los promedios
normal_promedio <- numeric(n)
gamma_promedio <-numeric(n)
poss_promedio <-numeric(n)

# Realizar el muestreo de las distribuciones.
for(i in 1:n) {
  normal <- rnorm(tamaño,mean = 5,sd =sigma/sqrt(n))
  normal_promedio[i] <- mean(normal)
  gamma <- rgamma(tamaño,alfa,rate=beta)
  gamma_promedio[i] <- mean(gamma)
  poss <- rpois(tamaño,lambda)
  poss_promedio[i] <- mean(poss)
}

# Código para gráficos
graficar_histograma <- function(medias_muestrales, titulo) {
  hist(medias_muestrales, probability = TRUE, main = titulo, xlab = "Media muestral")
  curve(dnorm(x, mean = mean(medias_muestrales), sd = sd(medias_muestrales)), col = "red", lwd = 2, add = TRUE)
}

graficar_histograma(normal,'Normal 1')
graficar_histograma(normal_promedio,'Normal con Promedio')
graficar_histograma(gamma,'Gamma 2')
graficar_histograma(gamma_promedio,'Gamma con Promedio')
graficar_histograma(poss,'Poisson 3')
graficar_histograma(poss_promedio,'Poisson con Promedio')
# Definición de variables o estructuras necesarias para el muestreo.
n <- 1000
tamaño <- 100
# Normal
sigma <- sqrt(5)
mu <-5
# Gamma
alfa <- 2
beta <- 1.5
# Poisson
lambda <- 5
normal_promedio <- numeric(n)
gamma_promedio <-numeric(n)
poss_promedio <-numeric(n)


# Realizar el muestreo de las distribuciones.
for(i in 1:n) {
  normal <- rnorm(tamaño,mean = 5,sd =sigma/sqrt(n))
  normal_promedio[i] <- mean(normal)
  gamma <- rgamma(tamaño,alfa,rate=beta)
  gamma_promedio[i] <- mean(gamma)
  poss <- rpois(tamaño,lambda)
  poss_promedio[i] <- mean(poss)
}
# Código para gráficos
graficar_histograma <- function(medias_muestrales, titulo) {
  hist(medias_muestrales, probability = TRUE, main = titulo, xlab = "Media muestral")
  curve(dnorm(x, mean = mean(medias_muestrales), sd = sd(medias_muestrales)), col = "red", lwd = 2, add = TRUE)
}
#Graficamos
graficar_histograma(normal,'Normal 1')

graficar_histograma(normal_promedio,'Normal con Promedio')

graficar_histograma(gamma,'Gamma 2')

graficar_histograma(gamma_promedio,'Gamma con Promedio')

graficar_histograma(poss,'Poisson 3')

graficar_histograma(poss_promedio,'Poisson con Promedio')

# Definición de variables o estructuras necesarias para el muestreo.
n <- 1000
tamaño <- 1000
# Normal
sigma <- sqrt(5)
mu <-5
# Gamma
alfa <- 2
beta <- 1.5
# Poisson
lambda <- 5
normal_promedio <- numeric(n)
gamma_promedio <-numeric(n)
poss_promedio <-numeric(n)


# Realizar el muestreo de las distribuciones.
for(i in 1:n) {
  normal <- rnorm(tamaño,mean = 5,sd =sigma/sqrt(n))
  normal_promedio[i] <- mean(normal)
  gamma <- rgamma(tamaño,alfa,rate=beta)
  gamma_promedio[i] <- mean(gamma)
  poss <- rpois(tamaño,lambda)
  poss_promedio[i] <- mean(poss)
}
# Código para gráficos
graficar_histograma <- function(medias_muestrales, titulo) {
  hist(medias_muestrales, probability = TRUE, main = titulo, xlab = "Media muestral")
  curve(dnorm(x, mean = mean(medias_muestrales), sd = sd(medias_muestrales)), col = "red", lwd = 2, add = TRUE)
}

#Graficamos
graficar_histograma(normal,'Normal 1')

graficar_histograma(normal_promedio,'Normal con Promedio')

graficar_histograma(gamma,'Gamma 2')

graficar_histograma(gamma_promedio,'Gamma con Promedio')

graficar_histograma(poss,'Poisson 3')

graficar_histograma(poss_promedio,'Poisson con Promedio')

# Definición de variables o estructuras necesarias para el muestreo.
n <- 1000
tamaño <- 5000
# Normal
sigma <- sqrt(5)
mu <-5
# Gamma
alfa <- 2
beta <- 1.5
# Poisson
lambda <- 5
normal_promedio <- numeric(n)
gamma_promedio <-numeric(n)
poss_promedio <-numeric(n)


# Realizar el muestreo de las distribuciones.
for(i in 1:n) {
  normal <- rnorm(tamaño,mean = 5,sd =sigma/sqrt(n))
  normal_promedio[i] <- mean(normal)
  gamma <- rgamma(tamaño,alfa,rate=beta)
  gamma_promedio[i] <- mean(gamma)
  poss <- rpois(tamaño,lambda)
  poss_promedio[i] <- mean(poss)
}
# Código para gráficos
graficar_histograma <- function(medias_muestrales, titulo) {
  hist(medias_muestrales, probability = TRUE, main = titulo, xlab = "Media muestral")
  curve(dnorm(x, mean = mean(medias_muestrales), sd = sd(medias_muestrales)), col = "red", lwd = 2, add = TRUE)
}
#Graficamos
graficar_histograma(normal,'Normal 1')

graficar_histograma(normal_promedio,'Normal con Promedio')

graficar_histograma(gamma,'Gamma 2')

graficar_histograma(gamma_promedio,'Gamma con Promedio')

graficar_histograma(poss,'Poisson 3')

graficar_histograma(poss_promedio,'Poisson con Promedio')

Pregunta 3: Ley de los Grandes Numeros.

Lanzamiento de monedas

Realice el experimento de lanzar una moneda cargada 1000 veces y observe el comportamiento que tiene la probabilidad de salir cara. Para realizar el experimento considere que la moneda tiene una probabilidad de \(5/8\) de salir cara. Grafique el experimento para las secuencias de intentos que van desde 1 a 1000, señalando el valor en que converge la probabilidad de salir cara.

Respuesta

convergencia <- 1:1000

# Simular lanzamientos
for (lanzamientos in 1:1000) {
  convergencia[lanzamientos] = rbinom(n = 1,size=lanzamientos,p=5/8)/lanzamientos
}

# Gráfico de la convergencia
#convergencia

plot(convergencia,type = "l", xlab = "Número de intentos", ylab = "Probabilidad de salir cara")
abline(h = 5/8, col = "red", lty = 2)

El problema de Monty Hall

Remontándonos en la televisión del año 1963, en USA existía un programa de concursos donde los participantes debían escoger entre 3 puertas para ganar un premio soñado. El problema del concurso era que solo detrás de 1 puerta estaba el premio mayor, mientras que detrás de las otras dos habían cabras como “premio”.

Una de las particularidades de este concurso, es que cuando el participante escogía una puerta, el animador del show abría una de las puertas que no fue escogida por el participante (Obviamente la puerta abierta por el animador no contenía el premio). Tras abrir la puerta, el animador consultaba al participante si su elección era definitiva, o si deseaba cambiar la puerta escogida por la otra puerta cerrada. Un vídeo que puede ayudar a comprender mejor este problema en el siguiente link.

Imagine que usted es participante del concurso y desea calcular la probabilidad de ganar el gran premio si cambia de puerta en el momento que el animador se lo ofrece. Utilizando listas/arrays/vectores simule las puertas del concurso, dejando aleatoriamente el premio en alguna posición del array. Hecho esto, genere un numero de forma aleatoria para escoger una de las puerta (posiciones de la estructura), para luego ver si cambiando de posición tendrá mayores posibilidades de ganar el premio. Genere N veces el experimento y grafique cada una de las iteraciones, tal como se hizo en el ejercicio de las monedas.

Respuesta:

# Creamos una función que simule el juego
montyhall <- function(cambiar = TRUE) {
  Puertas <- sample(1:3, 3) #Puertas donde la posición que tiene el 3 es el premio
  cat("Puertas:", Puertas, "; ")
  
  posicion <- sample(1:3, 1)  # Elección del participante
  cat("Posición elegida:", posicion, "; ")
  
  # Identificar puertas descartables que no tienen el premio y que no sean la elegida por el participante
  puertas_descartables <- which(Puertas != 3 & 1:3 != posicion)
  cat("Puertas descartables:", puertas_descartables, "; ")
  
  # Escoger una puerta para descartar
  puerta_descartada <- puertas_descartables[sample(length(puertas_descartables), 1)]
  cat("Puerta descartada:", puerta_descartada, "; ")
  
  # Si el participante decide cambiar su elección
  if (cambiar) {
    nuevas_opciones <- setdiff(1:3, c(posicion, puerta_descartada)) #se elige la otra puerta que queda disponible
    posicion <- nuevas_opciones
  }
  cat("Nueva elección:", posicion, "; ")
  
  # Verificar si la elección final tiene el premio
  Eleccion <- ifelse(Puertas[posicion] == 3, "Ganaste", "Perdiste")
  cat("Eleccion", Eleccion)
  return(Eleccion) #Retornamos la elección viendo si perdió o ganó
}

# Ejemplo de uso:
cat(montyhall(TRUE))


# Función que simula N juegos

n_juegos <- function(n = 10 ,cambiar_puerta = TRUE){
  convergencia <- numeric(n)
  
  for (juego in 1:n) {
    cat("Juego",juego , ":")
    resultado <- montyhall(cambiar_puerta)
    
    # Convertir el resultado a numérico: 1 para "Ganaste", 0 para "Perdiste"
    convergencia[juego] <- ifelse(resultado == "Ganaste", 1, 0)
  }
  return(convergencia)
}

# Ejecutar la simulación de 40 juegos cambiando de puerta
resultados_cambiar <- n_juegos(40, TRUE)

# Ejecutar la simulación de 40 juegos sin cambiar de puerta
resultados_no_cambiar <- n_juegos(40, FALSE)

# Graficar la convergencia para cambiar de puerta
plot(cumsum(resultados_cambiar) / 1:40, type = "l", col = "blue", xlab = "Número de juegos", ylab = "Proporción de éxitos", main = "Convergencia de Monty Hall", ylim = c(0, 1))

# Graficar la convergencia para no cambiar de puerta en el mismo gráfico
lines(cumsum(resultados_no_cambiar) / 1:40, col = "red")

# Añadir una leyenda para diferenciar las dos curvas
legend("bottomright", legend = c("Cambiar de puerta", "No cambiar de puerta"), col = c("blue", "red"), lty = 1)

Se puede ver como la elección de cambiar de puerta tiende a tener una mayor proporción de victorias que al no cambiarla.

Pregunta 4: ¿Independencia?

Ustedes disponen de los dados D1 y D2, los cuales no están cargados y son utilizados para comprobar que \(\mathbb{P}(AB)=\mathbb{P}(A)\mathbb{P}(B)\) cuando el evento A es independiente del B. Para estudiar la independencia considere que los eventos A y B se definen de la siguiente manera; sea A el evento dado por los valores obtenidos en el lanzamiento del dado D1, este está compuesto por \(A=\{D1=1,D1=2,D1=6\}\). Por otro lado, el evento B viene dado por los valores obtenidos con el dado D2, el que está conformado por \(B=\{D2=1,D2=2,D2=3,D2=4\}\). Con esto, tendremos un \(\mathbb{P}(A)=1/2\), \(\mathbb{P}(𝐵)=2/3\) y \(\mathbb{P}(AB)=1/3\). Compruebe de forma gráfica que al realizar 1000 lanzamientos (u otro valor grande que usted desea probar) se visualiza que \(\mathbb{P}(AB)=\mathbb{P}(A)\mathbb{P}(B)\).

Hecho lo anterior, compruebe el comportamiento de un segundo grupo de eventos, dados por el lanzamiento de solo el dado D1. Donde, los eventos para D1 quedan definidos como: \(A =\{D1=1,D1=2,D1=6\}\) y \(B=\{D1=1,D1=2,D1=3\}\). ¿Se observa independencia en este experimento?.

Se le aconseja que para simular los lanzamientos de dados utilice la función sample() para generar valores aleatorios entre 1 y 6. Compruebe los números generados por la función con los casos favorables de cada uno de los eventos a ser estudiados.

Primer grupo de eventos:

library(ggplot2)
# Primer grupo de eventos
N_lan = 1000 # Numero de lanzamientos

L_A = c(1, 2, 6) # Lanzamientos favorables A = c(1, 2, 6)
L_B = c(1, 2, 3, 4) # Lanzamientos favorables B = c(1, 2, 3, 4)
L_AB = c(1, 2) # Lanzamientos favorables AB = c(1, 2)

# Simular lanzamientos
D1 <- sample(1:6, N_lan, replace = TRUE)
D2 <- sample(1:6, N_lan, replace = TRUE)

A <- 0
B <- 0
AB <- 0

P_A_x_P_B <- numeric(N_lan)
P_AB <- numeric(N_lan)

# Conteo de casos
for (i in 1:N_lan) {
  a <- FALSE
  b <- FALSE
  if (D1[i] %in% L_A) {
    A <- A + 1
    a <- TRUE
  }
  if (D2[i] %in% L_B) {
    B <- B + 1
    b <- TRUE
  }
  if (a && b) {
    AB <- AB + 1
  }
  
  # Calcular las probabilidades acumuladas
  P_A_x_P_B[i] <- (A / i) * (B / i)
  P_AB[i] <- AB / i
}

# Gráfico del experimento
plot(1:N_lan, P_A_x_P_B, type = "l", col = "blue", ylim = c(0, 1), xlab = "Número de lanzamientos", ylab = "Probabilidad")
lines(1:N_lan, P_AB, col = "red")
legend("topright", legend = c("P(A)*P(B)", "P(A∩B)"), col = c("blue", "red"), lty = 1)

Respuesta:

Es claro observar que ambos eventos convergen a valores practicamente iguales a medida que aumentamos el número de lanzamientos , Esto confirma que el evento A y B son independientes, Ya que A depende del dado 1 y B depende del dado 2.

Segundo grupo de eventos:

# Segundo grupo de eventos
N_lan = 1000 # Numero de lanzamientos
  
L_A = c(1, 2, 6)  # Lanzamientos favorables A = c(1, 2, 6)
L_B = c(1, 2, 3)   # Lanzamientos favorables B = c(1, 2, 3)
L_AB = c(1, 2) # Lanzamientos favorables AB = c(1, 2)

# Simular lanzamientos
D1 <- sample(1:6, N_lan, replace = TRUE)

A <- 0
B <- 0
AB <- 0

P_A_x_P_B <- numeric(N_lan)
P_AB <- numeric(N_lan)

# Conteo de casos
for (i in 1:N_lan) {
  a <- FALSE
  b <- FALSE
  if (D1[i] %in% L_A) {
    A <- A + 1
    a <- TRUE
  }
  if (D1[i] %in% L_B) {  # D1 para B
    B <- B + 1
    b <- TRUE
  }
  if (a && b) {
    AB <- AB + 1
  }
  
  # Calcular las probabilidades acumuladas
  P_A_x_P_B[i] <- (A / i) * (B / i)
  P_AB[i] <- AB / i
}

plot(1:N_lan, P_A_x_P_B, type = "l", col = "blue", ylim = c(0, 1), xlab = "Número de lanzamientos", ylab = "Probabilidad", main = "Independencia - Un Dado")
lines(1:N_lan, P_AB, col = "red")
legend("topright", legend = c("P(A)*P(B)", "P(A∩B)"), col = c("blue", "red"), lty = 1)

Respuesta:

Es claro observar que ambos eventos no convergen a valores iguales , esto debido a que el evento A y B son dependientes, Ya que A depende del dado 1 y B tambien depende del dado 1.

Pregunta 5: La Ruina del Jugador

Un amigo ludópata suyo le comenta que el truco de jugar en el casino esta en no parar de apostar y apostando lo mínimo posible. Ya que así, tienes mas probabilidades de ganar el gran pozo que acumula el juego. Usted sabiendo la condición de su amigo, decide no creer en su conjetura y decide probar esto a través de un experimento.

Para realizar el experimento usted decide asumir las siguientes declaraciones, bajo sus observaciones:

En el primer experimento deberá obtener la evolución de los fondos hasta que el jugador se queda sin fondos para jugar. Puede ser útil seguir la lógica de una moneda cargada para realizar esto. Pruebe esto con una apuesta igual a 5, 25 y 50 graficando los resultados. Comente los resultados obtenidos.

Para la segunda parte de este experimento, con las funciones ya creadas, proyecte 5000 apuestas y obtenga la probabilidad a la que converge el experimento empezando con una apuesta de: 5, 25 y 50. Para este experimento considerar como éxito (1) cuando su función ruina supera los 200 y considere lo contrario cuando su función perdida es menor o igual a 0.

Respuesta

# Función para obtener el desarrollo de las apuestas
ruina <- function(fondos = 100, apuesta = 5){
  vec_fondos <- c(fondos)  # Monto inicial
  prob_ganar <- 8 / 19    # Probabilidad de ganar
  while (0<fondos & fondos<=200) {
    if (runif(1) < prob_ganar) {  # Generar un número aleatorio para determinar el resultado del juego
      fondos <- fondos + apuesta  # Ganó
    } else {
      fondos <- fondos - apuesta  # Perdió, disminuye los fondos
    }
    vec_fondos <- c(vec_fondos, fondos)  # Guardar los fondos despues de cada iteracion
  }
  return(vec_fondos) # Devuelve un vector con el desarrollo de los fondos
}

plot(ruina(), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 5)")
plot(ruina(apuesta = 25), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 25)")
plot(ruina(apuesta = 50), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 50)")

Respuesta :

Con una apuesta de $5, la evolución de los fondos muestra una tendencia a la baja, el jugador juega muchas veces antes de quedarse sin dinero. Es difícil que el jugador supere los $120, ya que las ganancias pequeñas no compensan las pérdidas acumuladas.

Con una apuesta de $25, el jugador se queda sin fondos en menos tiempo comparado con una apuesta de $5.

Con una apuesta de $50, las pérdidas se sienten de manera aún más rápida, pero las ganancias, cuando ocurren, son mayores, llegando a ganar considerablemente más veces comparado con los otros dos casos.

ruina <- function(fondos = 100, apuesta = 5){
  vec_fondos <- c(fondos)  # Monto inicial
  prob_ganar <- 8 / 19    # Probabilidad de ganar
  while (0<fondos & fondos<=200) {
    if (runif(1) < prob_ganar) {  # Generar un número aleatorio para determinar el resultado del juego
      fondos <- fondos + apuesta  # Ganó
    } else {
      fondos <- fondos - apuesta  # Perdió, disminuye los fondos
    }
    vec_fondos <- c(vec_fondos, fondos)  # Guardar los fondos despues de cada iteracion
  }
  return(vec_fondos) # Devuelve un vector con el desarrollo de los fondos
}


probabilidad_ruina <- function(apuesta = 5, n_simulaciones = 5000){
  exitos <- 0  # Contador de éxitos
  
  for (i in 1:n_simulaciones) {
    final_fondos <- ruina(apuesta = apuesta)
    if (final_fondos[length(final_fondos)] > 200) {
      exitos <- exitos + 1  # Exito si se superan los 200 dólares
    }
  }
  
  prob_exito <- exitos / n_simulaciones  # Calcular la probabilidad de éxito
  return(prob_exito)
}

# Calcular la probabilidad de ruina para diferentes tamaños de apuesta
prob_ruina_5 <- probabilidad_ruina(apuesta = 5)
prob_ruina_25 <- probabilidad_ruina(apuesta = 25)
prob_ruina_50 <- probabilidad_ruina(apuesta = 50)

cat("Probabilidad de ruina con apuesta de 5 dólares:", prob_ruina_5, "\n")
## Probabilidad de ruina con apuesta de 5 dólares: 0.001
cat("Probabilidad de ruina con apuesta de 25 dólares:", prob_ruina_25, "\n")
## Probabilidad de ruina con apuesta de 25 dólares: 0.1548
cat("Probabilidad de ruina con apuesta de 50 dólares:", prob_ruina_50, "\n")
## Probabilidad de ruina con apuesta de 50 dólares: 0.2182

A work by CC6104